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Abstract 

This paper explores solutions to the spherically symmetric Euler equations. 
Motivated by the work of Hagstrom and Hariharan [7] and Geer and Pope [5], 
we model the effect of a pulsating sphere in a compressible media. The literature 
available on this suggests that an accurate numerical solution requires artificial 
boundary conditions which simulate the propagation of nonlinear waves in open 
domains. Until recently, the boundary conditions available are in general linear, 
and based on non-reflection. Exceptions to this are the nonlinear non-reflective 
conditions of Thompson [11], and the nonlinear reflective condition of [7]. The 
former is based on the rate of change of the incoming characteristics, while the 
latter relies on asymptotic analysis and the method of characteristics and accounts 
for the coupling of incoming and outgoing characteristics. Furthermore in [7] it 
was shown in a test situation in which the flow would reach a steady state over 
a long time, the method proposed in [11] could lead to an incorrect steady state. 
The current study considers periodic flows. Moreover, all possible types and 
techniques of boundary conditions are included in this study. The technique 
recommended by [7] proved superior to all others considered, and matched the 
results of asymptotic methods which are valid for low subsonic Mach numbers. 

"Supported in part by NSF Grant DMS-8921189, and by NASA cooperative 
agreement NCC-3-104. 



1 Introduction 


This paper deals with an exterior problem governed by compressible, in viscid gas dynam- 
ics equations. The nonlinear nature of the problem calls for a numerical solution to the 
problem. As common to most exterior problems, the condition(s) that must be satisfied 
at infinity is the key to solve these problems uniquely. The conditions that are physically 
correct are not necessarily mathematically correct for the well-posed-ness of the problem. 
Among the well known examples the one that stands out in the literature is the exterior 
problem governed by the Helmholtz equation. The class of problems that is governed by 
this equation include scattering of acoustic and electromagnetic waves. While physically, 
the quantities decay to zero at infinity, the requirement being exactly zero at infinity is 
not a well posed condition, as observed by Hefiwig (see referenceflO]). The same holds 
true for the Euler equations under consideration. Moreover, the computational consider- 
ations require that these problems are to be solved in a finite domain with the aid of an 
artificial boundary placed sufficiently far away Horn the region of local flows that are of 
interest. This introduces another parameter L (eg. the radius of the artificial boundary) 
and the requirement of boundary condition(s) that should be consistent with the true 
infinite domain problem. Several efforts have been made along these lines. Most efforts 
consider linearized flows (linearized about the state at infinity) to obtain such classes of 
boundary conditions (For example see [2] and [6]). Among the first nonlinear versions of 
the boundary conditions for gas dynamics problems was that of Hedstrom [9]. Recently 
Thompson [11] extended the ideas of Hedstrom to multidimensional problems. His pro- 
cedure works well for two dimensional problems described in Cartesian coordinates with 
rectangular boundaries. Thompsons method works because in Cartesian coordinates the 
flux term of the equation for the outgoing wave can be neglected. Although this technique 
can be implemented in spherical and cylindrical coordinate systems, it will produce poor 
results because the equations for the outward and inward propagating wave variables are 
coupled in the flux term. 

There is another parameter that is present in the computational model. This param- 
eter is the total length of computational time, T. For problems that require long time 
stability, such as the ones considered in this paper, the ideas of [3] cannot be used as 
these are high frequency boundary conditions. They will be accurate for short times and 
one cannot expect them to work well for any long time calculations. 

Motivated by the work of Hagstrom and Hariharan [7] and the asymptotic work of 
Geer and Pope [5] as well as Hardin and Pope [8], we model the external flow around a 
sphere pulsating in the air. This problem is of considerable interest from a mathematical 
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standpoint as a testing ground for various versions of time dependant open boundary 
conditions and asymptotic methods. Until recently, such boundary conditions have all 
been linear, and based on non-reflection, however for many problems such as the pulsating 
sphere, or jet flow computations, the the nonlinearity is important, and the waves are 
coupled such that an outward propagating wave will generate an inward travelling wave. 
Exceptions to these kind of conditions are the technique of [7], which is a nonlinear 
reflective condition that is based on asymptotic analysis, and the nonlinear non reflective 
technique of Thompson [11]. This work applied the above methods to the case of fully 
unsteady flow resulting from a pulsating sphere. Furthermore various other boundary 
conditions were assessed along with the ones in [7] and [11]. 
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2 Model and Derivation 


Considering isentropic, compressible flow the Euler equations are given by 

^ + V • (pS) = 0, 


( 1 ) 


A»i 

Vp + p(^ + (ff' v )i)=°, (2) 

V = kp y . (3) 

Equations (1) and (2) are known as the conservation of mass and momentum equations, 
and (3) is an isentropic constitutive law. Here p is density, u is the velocity field vector, 
p is pressure, t is time, and k and 7 are constants characteristic to the medium. The 
model problem considered is the special case of a sphere harmonically pulsating in a fluid 


medium. Accordingly, the first assumption in this problem is 
symmetric and the governing equations take the form 

that the flow is spherically 

dp d . . 2pu 

(4) 

d , . d ( 2 \ 2pu 2 

at (r*)+ gr (p + im)= r , 

(5) 

p = kp 1 . 

( 6 ) 

Then we scale the variables in the following manner: 


II 

$ 

35 . 

(7) 

u = Uu , 

(8) 

r = Irf, 

(9) 

t = ( L/Co)t , 

( 10 ) 


where U, po, and L are a radial velocity, density, and length characteristic to the problem 
and the local speed of sound of the medium is defined by 


c* = ^ (11) 

dp 

Co = c(po), (12) 

with the Mach number given as M — U/co. This scaling along with letting z — Mpu 
and implementing ( 6 ) yield the following nondimensionalized system 


dp dz —2 z 

dt dr r ’ 


(13) 
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+ + = (14) 

dt dr \ 7 p ) rp 

Desiring that the sphere pulsate from a radius of unity at time t = 0 the pulsation of 
frequency w is implemented as a boundary condition on it’s surface, 

tt(l ,<) = Main(wt) } (15) 

while keeping the surface of the sphere at r = 1 to simplify computational aspects of the 
problem. Density, which is not specified on the sphere’s surface, will be recovered from 
the governing equations through the method of characteristics. In the far field the gas 
should be undisturbed, thus p -► 1, and v -> 0 as r -► oo. Initially we want the flow at 
rest, which implies that the initial conditions are p(r, 0) = 1 and u(r, 0) = 0. 
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3 Asymptotic Analysis 


It was shown in [7] that the Euler equations (13) and (14) can be cast in terms of Riemann 
variables, R and S t as follows: 


where 


dR i 


dR _ 

2 i? 7-3 

*i 

+ 

:p + p ’) 

dr 

~T P 7 

as , 

( z 

1 as 

2 2 2=1 



' dr ~ 

T f ’ ’ 


R = - + 

2 p 3 ? 



P 

7-1' 



(16) 

(IT) 

(18) 


and 


P 

This form is essential to examine the structure of far field solutions. Recall that the 
boundary conditions of the pulsating sphere axe: u(l,t) = Msin(wt), and also p — * 1, 
and if — > 0 as r — ► oo. Unfortunately an infinite domain is impossible to model in a 
computational sense, thus artificial boundary conditions are demanded. One method 
is to approximate infinity by a very large number, say L, and prescribe p(L,t) = 1 
and tt(£, t) = 0. Although simple to implement, this approximately physically correct 
boundary condition is not recommended since the effect of the pulsation must never 
propagate all the way to L or the solution will become invalid since the state variables are 
still time dependant for any finite L. The preferred way of handling an infinite domain is 
to adopt an Eulerian framework and observe the flow as it passes through a fixed volume. 
There are several approaches to building such boundaries. The approach recommended 
in this discussion will stem from an asymptotic analysis of the Riemann variables and will 
closely mimi c the technique of [7] which proposed that if R and S could be approximated 
on the boundary then a more accurate artificial boundary condition could be developed. 
Thus we expanded the wave variables radially. 


2 pV 
7-1 


(19) 


R = Ro + 


-Ri(M) 

r 


Rifat) 
^ *.2 


+ ••• 


( 20 ) 


s = So + ^) + ^i) + 


OO 


r r* 

Using the fact that R, and S are given in terms of the fluid properties 


R = \~ Gy 

p 


( 21 ) 

( 22 ) 
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(23) 


S = - - G, 

P 

the values of Rq, and S 0 come from direct application of the far field data. 



C* 1 

II 

to 

1 

II 

(24) 

Then by substituting 

^ = (V) <* - s) ' 

(25) 

and 

z R + S 

P ~ 2 1 

(26) 

the wave variable equations assume the following form, 


“ + ( 
dt ' 

'R\S 7-lR-S\8«_ 

, 1 11 S 1 2 J Jr » ' ' 

(27) 

" + | 
dt ' 

rJJ+S 7-1B-SN3S 

<2 1 2 j 2 ) dr 2 r V > 

(28) 


Applying the asymptotic expansions 

^(Rjr + Rt/r* + •••) + 


(Ri + Si)/2r + ( R 2 + S 2 )/2r 2 + • 

+(7 - l )/2 {Ro + (Ri - 5i)/2r + (R 2 + 5 2 )/2r a + •■•)] 

= {2{R 0 R l - S 0 S 1 )/r + (2 RqR 2 + R\ - 2 S 0 S 2 - S*)/r 2 + •••), 


(Ri/t + • • •) (29) 


and, 


£(Si/r + S,/r 2 + •••) + 
(R l + S\)/2r + ( R 2 + S 2 )/2r 2 + 


^ (*1 + Si)/2r + (R 2 + S 2 )/2r + • • • JL ( Sl fr + . . •) 

+ (7 - l )/2 (Ro + (Ri ~ 5x)/2r + (R 2 + S 2 )/2t 2 + • • *) J * V 

= -atl (2(J2 oJ2i - S 0 Si)/r + (27^2 + Rl~ 2 S 0 S 2 - S 2 )/r 2 + •••)> 


we obtain the first order problem, 

dR\ ( R\ + 7 ~~ j- 


R\ + Si \ \ dR\ 


■))$-* 


m fR±±Si + iz±{ R 

dt \ 2r 2 V 

dS 2 , (Ri + Sr 7 - 1 (j, . Ri + S A\ - n 

-ar + H; r *• - 0 ’ 


(30) 


(31) 

(32) 
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which fits the form of Riemann invariants, that is, along the proper characteristic curves, 
Ri and Si are constants. Furthermore, since there should be no incoming waves at t = 0 
then S\(r, t) = 0. The characteristic curve for R\ given by 


dr 

H 



Ri 

r 


(33) 


defines the slope of the paths through space- time on which Ri is constant, so that Ri can 
be written as a function of only one variable distinguishing which curve it came from, 
H(t) = /Ji(r,i(r,r)). Thus 


/ 


dr 


1 + (i + *?) 



(34) 


Then requiring that t = r when r = L, 

= t — r. 

Note that we have not yet calculated R\ on the boundary. The proposed boundary 
condition actually follows manipulation of the 5 j equation. 



L - + pp) tf(*01og 


(j + 3 T i ) ^Hi- 


fi + 3 7 1 )ff( 1 -)+ i 


dSi 

dt 



(36) 


Writing this in in terms of the characteristic for Ri 




= H(r) 

(37) 

as 2 

dr 



(38) 

But 

?L = i 
dt 


(39) 

and 

dr ( 7 + 1 )J5T(r) 

dr 4r 2 + (7 - 3)tf(r) 


(40) 
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so that 


Now let 


and 


it [ J " (* + tL ir^^ (r) ) l 1 “ 4r>+h - 3)tf(r))] 
X,N(r,r) = 1 - (l + ^*) (l - 4r ,i 7 (7 + _ 1) 3 ^ (r -j) 


D(r,r) = ^1 + ^ 7 4r ~ g ( r )^ • 
Equation (38) can then be rewritten as 


^ - N(t, t)D{t,t)^~ = N(t,t)H(t). 


Now along dr I dr = — iV(r,r)D(r,r) 


dr 

Given the initial conditions 52—0 for t = 0, 52 at L is given by 

S,(L,t)= £ H(.)N(>,((r,T,L))d. 

Jo 

where 


d£ 

ds 


= -ND, 


and 


t(r;T,L) = L. 

Hence upon differentiating with respect to r 


dSj 

dr 


= H(t)N(t , L) + £ H(s)^N( 3, L))ds 


~ = Ri(L,t)N(t, L) + £ H(.)£ ( N(.,((.‘-,r,L))^i>. 


(41) 

(42) 

(43) 

(44) 

(45) 

(46) 

(47) 

(48) 

(49) 

(50) 
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Then considering 


8N(t, r) 

dr 


(7 + l)8rg(r) 2(7 - 1)<7 ~ 3)g» , 7 ~ 3 , 

(4r> + (7 - 3)tf(r)) ! (4r> + (7 - 3)g(r)) a 4r> 


(7 + 0(7 ~ W) 

4r a (4r* + (7 - 3 )H(r)) 


X 




(7+l)g 

4i* 3 + (7 — 3)g(r) / / ' 


( 51 ) 


which implies that dN/d£ causes the integral in (50) to be of order 1/r and thus con- 
tributing to the terms S3 and higher. Using R = Ro + & = ^0 + ^*^2 the 

boundary condition at r = L becomes, 


es m.t )-Jj 

8t . 2 L 

T—Lt 


In terms of the primitive variables we have 

2L 


(52) 


(53) 


Thus the physical significance of this is that in the far field there is no incoming wave, 
and that at any fixed distance L the correction will be given in terms of the outgoing 
wave. 
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4 Numerical Implementation 


4.1 Numerical Scheme 

We have only considered a relatively simple second order accurate Lax-Wendroff algo- 
rithm. Higher order schemes can be applied with care in handling shock waves that are 
present in the current problem. To begin with, the continuity and momentum equations 
have been presented in a weak conservation law form, U t + F r = — W, where, 


U = 



(54) 



In this format any of several schemes discussed by G.A. Sod [12] can be implemented. 
The notational convention adopted for the purpose of discretization is that a function of 
space- time, F (r,t) = F”, with radius discretized by *, and time by n. In keeping with 
[7], the two step Lax-Wendroff scheme was chosen for its smoothness, and simplicity. For 
programming purposes the U vector is written as 



The half step time values of U are denoted by V“ . The updated values of UJ* can then 
be written right back on top of U" to save memory. The first step is given by 


v? = i (u? + u? +1 ) - 5& (F(U” +1 ) - F(U" )) - ¥ W(1(U? + U? +1 W}> 

; = !,•••, n 


(58) 


and the second step is, 

u? +1 = U? - ^ (FfV?) - F(V?.,)) - A(W (i(VT - Vfc,),n) (59) 

The next step is to implement the flux corrective transport algorithm (FCT) of Fletcher 
[4] in order to analyze the data and remove any spurious oscillations generated by the 
presence of shock waves. 
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4.2 Initial and Boundary Conditions 

Implementation of this scheme requires initial and boundary data. Initially U(l, t, 1) = 1, 
and U(2,», 1) = 0. The boundary data will come from the boundary conditions chosen 
at r = 1 and r = L which are discretized as * = 0 and i = N + 1 respectively. We will 
discuss five different methods. 

4.2.1 Method One 

Method one is the physically approximate conditions taken from the ambient state in the 
far field. 

p(L,t) = 1 

and 

z(I,t) = 0 


4.2.2 Method Two 

Method two is Thompson’s technique [11], which is to neglect the flux term in the gov- 
erning S equation. This means at r = L we have 


dS _ 2 zp 1 ^ 
~dt = r ' 


(60) 


In discretized form it reads 

STO+x) = sm + 1) + S(tft) - S(Uff +1 ) + 2A tzN+rps+^/L. 

Then use a discretized version of (71) to update R. 


4.2.3 Method Three 

Method three totally suppresses the incoming wave variable, that is set 



Discretized for the nonlinear case, update S as 

S(tO = S(Eft +1 ) + S(t7S) - S(Uj?') 

then, xis in method two, use (71) to update R. 


(61) 


(62) 


(63) 
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4.2.4 Method Four 

Method four is the linearized version of method three, that is take 



R = — 2 — + (p-l) + *, 

7-1 

(64) 


S = — r - (p - 1) + ^, 

7-1 

(65) 

and then apply 

SIS 

n 

o 

(66) 


and equation (71) to update R and S. 


4.2.5 Method Five 


The boundary condition which we will show to have the best results is the one proposed 
in [7]. 

as 

dt 


r=L 


2L 


(67) 


where L is the truncation radius. Use of this boundary condition calls for an updated 
value of R. This means that not only (67) must be discretized, but also the partial 
differential equation for R. Writing (67) in the form 


as 

dt 


Q(p, z ) 


( 68 ) 


A simple second order discretization is, 

S(Un\\) = S{Ujf+i) + S(V%) - S{UZ +1 ) + 2AtQ(V£) 

where i = N -f 1 corresponds to r = L. Next the R equation in the form 


dR s dR . 

-^ + C{p,z)-^-H(p,z) 


(69) 


(70) 


is discretized to, 


RMp = 

(i/(i + gc(vs)) [fl(£/» +1 ) + JJTC) 
-R(Utf') + 2Atff(V?‘) 

-£,CW) (W«) - RW) - W))] 


( 71 ) 
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Then the primitive variables can be updated at the boundary using the updated boundary 
values of R and S with , 


(V ( * s) )’ ■ 

( 72 ) 

z = 2 ^ + 

( 73 ) 
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5 Results and Comparisons 

As mentioned earlier, the governing equations were discretized using the two step Lax- 
Wendroff scheme. All numerical tests were conducted with the surface of the sphere 
pulsating with radial velocity u = Msin(wt) where w = 1.5, and the Mach number set 
at, M = 0.5. The tests were designed to produce radial density and momentum profiles 
at time value t = 10. Initially we ran tests using a computational domain that was larger 
than the radius of propagation. This guaranteed that we had a solution that was inde- 
pendent of any boundary conditions. With a Mach number of M = 0.5 Gibbs phenomena 
became visible near the shock waves, displaying highly oscillatory behavior. A flux cor- 
rective transport subroutine by Fletcher [4] was employed to help capture the shocks and 
minimize oscillations as they occurred. To reassure the validity of the numerical solution 
it was necessary to compare it with a completely independent method. Asymptotics are 
often considered as alternatives to numerical solutions. Asymptotic solutions use the 
Mach as a parameter of expansion, thus requiring very low Mach numbers. The multiple 
scales technique of Geer and Pope [5] was chosen to compare with method one. Their 
first order approximation for velocity is given as 


«o(r,t) 


(r — l)tt7 
(1 + tu 2 )r s 


cos(r) + 


1 

(l+u» a ) 


(w 2 /r + 1/r 2 ) sin(r) 


and the second order correction is given as 


where 


and 


«i(r, t) = (&a/r J + 26at o/r) cos(2r) + ( h/r 2 - 2b2w/r) sin(2r) 


( (29 - 37)u> 2 w 2 {w 2 - 1)\ . (l+7)u; 2 log(r) 

\16r 3 (l + w 2 ) 2 + r 2 (l + w 2 ) 2 ) C °^ ^ 4r 2 (l + w 2 ) 


( (29 - 37)w(u> 2 — 1) 
^ 32r 3 (l + w 2 ) 2 


2w 3 

r 2 (l + tn 2 ) 2 



f(t) = (sin(urf) - w cos(t vt)/ (1 + w 2 ) 
r 2 = r, 2 + (7 + 1 )M{t - 1 )f\t + 1 - v) 

T = w(t + 1 — 7 ) 


3u> 2 (14 - 27 + 17u> 2 + 7 w 2 ) 

! = ” 16(1 + tn 2 ) 2 (l + 4w 2 ) 

w(29 — 37 — 17w 2 + 157tn 2 — 64u/ 4 ) 
32(1 + u> 2 ) 2 (l + 4to 2 ) 


(74) 


(75) 


(76) 

(77) 

(78) 

(79) 

(80) 
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Then the velocity is calculated as 

u(r, t) = vo(t } t) + M«i(r, t). (81) 

Figure 1 shows a velocity profile of [5] compared with the numerical method at time, 
t = 10, generated with Mach number, M - 0.05. These solutions agree fairly closely. 

If an artificial boundary was placed at a radius where waves would pass through it 
by time t = 10 then the effect of the artificial boundary conditions could be evaluated 
by comparing it to the profiles calculated in a large domain which was independent of 
boundary conditions. A perfect boundary condition would yield a profile that perfectly 
matched these examples regardless of where the radius was truncated. In figure 2 the 
physically approximate boundary conditions of method one were evaluated at r = 6, 
r = 7, and r = 9. Figure 2 clearly shows that this physical approximation is invalid for 
truncations inside of the radius of propagation. 

In figure 3 Thompsons method [11], method two, has been run at three different 
truncation radii, L = 21 (X = 21 corresponds to 2000 nodes.), L — 7, and L = 3. If the 
method two boundary condition was numerically correct, the profiles for L = 7, and L = 3 
should virtually coincide with the profile for L = 21. Thompson’s boundary condition 
Ting caused the numerical scheme to undershoot the L = 21 profile when truncated at 
L = 3, this is because Thompsons boundary condition will go to zero in the far field 
independent of the behavior of the density. 

The results of method three, the nonlinear suppression of the incoming wave, also 
truncated at L = 21, L = 7, and L = 3, are given in figure 4. The basic difficulty is 
the fact that even though this method was based on the Biemann variables, it did not 
simulate the movement of the waves because it did not exploit the governing equations 
of the Riemann variables (16) and (17). 

Figure 5 presents the results of method four, the linearized suppression of the incoming 
wave, when truncated at L = 21, L = 7, and L = 3 as in the previous tests. Although 
not perfect, this technique does not perform badly for ten time units. The amplitude of 
overshoot from the L = 3 truncation indicates that this method will not stay stable in a 
far time application. Thus as time is increased the profile would continue to deviate. 

The results of the technique proposed in [7], method five, as shown in figures 6, 7, 
and 8, are far superior to the other four methods which offer no hope for solutions in 
long time calculations . In figure 6 the profiles generated with truncations at L = 3 and 
L = 7 lay precisely on the profile of truncation at L = 21. This same style test was then 
performed on method five at time t = 100 and t = 1000 (one million time steps) to see if 
it would re main stable for a far time situation. In figures 7 and 8 the truncated profiles 
continue to coincide, indicating that this technique precisely simulates waves passing the 
artificial boundary and is dependable regardless of the time value to which the scheme is 
run. 
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It is useful to examine the nature of the solutions as w -> 0 and as w -*■ oo. As 
one would expect, the variation of the solution in the computational domain are small 
for low frequencies. Figure 9 shows this for decreasing values of w using method five. 
(In this case the nonlinear effects are more prominent in the far field than in the near 
field.) On the other hand for high frequency waves an increase in the frequency reflects 
more periodic structures of the shock waves. Figure 10 shows a tenfold increase in the 
number of shock fronts by going from w — 1.5 to w = 15. However further increase 
in the frequency poses a challenge for the numerical computations. As the frequency 
increases the fronts are closer to each other. In the extreme cases one must develop an 
asymptotic solution to these problems. For example figure 10 suggests that the solution 
rapidly attenuates for w = 150. These asymptotic aspects require further study. 

Considering figure 9, method five works quite well for low frequency problems, but 
this was not the case for all other classes of boundary conditions. In figure 11 we show 
the solution of Thompson for w — .015, with L — 11 and t = 100. For t = 100 the 
sphere will not have gone through a full oscillation, but the wave will have passed our 
truncation point at L = 11. This solution approaches an incorrect state well below what 
would be expected (compare with figure 9), Thus it is not uniformly accurate in time. 
The method proposed in [7] is indeed a uniformly valid condition due to its derivation 
that is based on uniform asymptotic expansions. 
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6 Conclusion 


This paper has modeled the flow field generated by a pulsating sphere. This model- 
ing process involved exploring the laws of gas dynamics and combining these laws with 
mathematical analysis in order to obtain numerical solutions for the flow field. Because 
of the intrinsic spherical symmetry of the problem the solutions naturally took the form 
of radial density and velocity profiles for specified times of interest. The numerical im- 
plementation of the governing equations was of particular mathematical interest because 
it became necessary to include a flux corrective transport algorithm to help control the 
numerical instabilities generated by the presence of shock waves. 

The next topic to be addressed in the modeling process was to come up with suitable 
numerical boundary conditions. Following the lead of Hagstrom and Hariharan [7] a 
method of characteristics representation was created to sense the presence of incoming 
and outgoing waves, and then implemented as a numerical boundary condition. The 
boundary condition was then evaluated along with several other techniques for stability 
far in time, and consistency regardless of the radius at which the domain was truncated. 
This discussion revealed that for the numerical solution of problems in gas dynamics the 
modeling of the open boundary conditions is critical. The wrong form of these conditions, 
especially in fully unsteady flow, can lead to a completely incorrect solution. Yet from 
all this, the technique of [7] has proven to be a dependable method for problems that 
can be modeled as having spherical wavefronts in the far field. Future work in this field 
is the generalized boundary condition for fully three dimensional unsteady flow. 
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Figure 1: Asymptotic technique vs. Numerical method for Mach number M = 0.05 



Figure 2: Physically approximate boundary conditions 
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Figure 3: The boundary condition of Thompson 



Figure 4: Nonlinear suppression of the incoming wave 
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Figure 6: The method of Eagstrom and Hariharan 
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Figure 9: Density profiles from low frequency pulsation, t=100 



Figure 10: Density profiles from high frequency pulsation, t=100 
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Figure 11: Density profiles for w— .015, r— 11, and t— 100 
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